{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>610</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>613</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>613</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     z\n",
       "0  610\n",
       "1  613\n",
       "2  613"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>640</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>640</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>640</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       z\n",
       "996  640\n",
       "997  640\n",
       "998  640"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd   ###variationalinference1init\n",
    "import math, random\n",
    "all_data  = pd.read_csv(\"sensor_data_600.txt\", delimiter=\" \", header=None, names = (\"date\",\"time\",\"ir\",\"z\"))#lidarのセンサ値は「z」に\n",
    "data = all_data.sample(1000).sort_values(by=\"z\").reset_index()  #1000個だけサンプリングしてインデックスを振り直す\n",
    "data = pd.DataFrame(data[\"z\"])\n",
    "\n",
    "display(data[0:3], data[-4:-1]) #とりあえず最初と最後のデータを表示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>610</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>613</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>613</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     z    0    1\n",
       "0  610  1.0  0.0\n",
       "1  613  1.0  0.0\n",
       "2  613  1.0  0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>640</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>640</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>640</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       z    0    1\n",
       "996  640  0.0  1.0\n",
       "997  640  0.0  1.0\n",
       "998  640  0.0  1.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "##負担率の初期化## ###variationalinference1rate\n",
    "\n",
    "K = 2 #クラスタ数\n",
    "n = int(math.ceil(len(data)/K)) #クラスタあたりのセンサ値の数\n",
    "for k in range(K):\n",
    "    data[k] = [1.0 if k == int(i/n) else 0.0 for i,d in data.iterrows()] #データをK個に分けて、一つのr_{i,k}を1に。他を0に。\n",
    "    \n",
    "display(data[0:3], data[-4:-1])  #下の出力の「0」、「1」が負担率の列"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_parameters(ds, k, mu_avg=600, zeta=1, alpha=1, beta=1, tau=1):    ###variationalinference1params\n",
    "    R = sum([d[k] for _, d in ds.iterrows()])\n",
    "    S = sum([d[k]*d[\"z\"] for _, d in ds.iterrows()])\n",
    "    T = sum([d[k]*(d[\"z\"]**2) for _, d in ds.iterrows()])\n",
    "    \n",
    "    hat = {}\n",
    "\n",
    "    hat[\"tau\"] = R + tau\n",
    "    hat[\"zeta\"] = R + zeta\n",
    "    hat[\"mu_avg\"] = (S + zeta*mu_avg)/hat[\"zeta\"]\n",
    "    hat[\"alpha\"] = R/2 + alpha\n",
    "    hat[\"beta\"] = (T + zeta*(mu_avg**2) - hat[\"zeta\"]*(hat[\"mu_avg\"]**2))/2 + beta\n",
    "    \n",
    "    hat[\"z_std\"] = math.sqrt(hat[\"beta\"]/hat[\"alpha\"])\n",
    "    \n",
    "    return pd.DataFrame(hat, index=[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>tau</th>\n",
       "      <th>zeta</th>\n",
       "      <th>mu_avg</th>\n",
       "      <th>alpha</th>\n",
       "      <th>beta</th>\n",
       "      <th>z_std</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>501.0</td>\n",
       "      <td>501.0</td>\n",
       "      <td>622.093812</td>\n",
       "      <td>251.0</td>\n",
       "      <td>3046.295409</td>\n",
       "      <td>3.483767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>501.0</td>\n",
       "      <td>501.0</td>\n",
       "      <td>631.918164</td>\n",
       "      <td>251.0</td>\n",
       "      <td>2604.822355</td>\n",
       "      <td>3.221456</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     tau   zeta      mu_avg  alpha         beta     z_std\n",
       "0  501.0  501.0  622.093812  251.0  3046.295409  3.483767\n",
       "1  501.0  501.0  631.918164  251.0  2604.822355  3.221456"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = pd.concat([update_parameters(data, k) for k in range(K)])   ###variationalinference1paramsolve\n",
    "params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, dirichlet ###variationalinference1draw\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def draw(ps):\n",
    "    pi = dirichlet([ps[\"tau\"][k] for k in range(K)]).rvs()[0]\n",
    "    pdfs = [ norm(loc=ps[\"mu_avg\"][k], scale=ps[\"z_std\"][k]) for k in range(K) ]\n",
    "\n",
    "    xs = np.arange(600,650,0.5)\n",
    "\n",
    "    ##p(z)の描画##\n",
    "    ys = [ sum([pdfs[k].pdf(x)*pi[k] for k in range(K)])*len(data) for x in xs] #pdfを足してデータ数をかける\n",
    "    plt.plot(xs, ys, color=\"red\")\n",
    "\n",
    "    ##各ガウス分布の描画##\n",
    "    for k in range(K):\n",
    "        ys = [pdfs[k].pdf(x)*pi[k]*len(data) for x in xs]\n",
    "        plt.plot(xs, ys, color=\"blue\")\n",
    "\n",
    "    ##元のデータのヒストグラムの描画##\n",
    "    data[\"z\"].hist(bins = max(data[\"z\"]) - min(data[\"z\"]), align='left', alpha=0.4, color=\"gray\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4lFXa+PHvmfQeUhgCREB6DSSht4A0UQR7W2XXwuoWddf6vrrvttd9d/25Kuuuu8vaWEVRcVFExIKELpAECL0TCJCEBEJ6mcz5/TETDJCQyWT63J/rypWZZ55yn+SZe545zylKa40QQgjvZ3B3AEIIIRxDEroQQvgISehCCOEjJKELIYSPkIQuhBA+QhK6EEL4CEnoQgjhIyShCyGEj5CELoQQPiLQlQdLSEjQ3bt3t2vbyspKIiIiHBuQF/DHcvtjmcE/yy1ltk12dnax1jqxtfVcmtC7d+9OVlaWXdtmZmaSkZHh2IC8gD+W2x/LDP5ZbimzbZRSebasJ1UuQgjhIyShCyGEj5CELoQQPkISuhBC+AhJ6EII4SMkoQshhI+QhC6EED5CEroQQvgISehCCOEjXNpTVAhhm+zs7GaXp6WluTgS4U1sukJXSv1CKbVbKbVLKfW+UipUKdVDKbVZKXVIKfWBUirY2cEKIYRoWasJXSnVBXgESNdaDwICgDuAPwEva617AeeA+50ZqBBCiCuztQ49EAhTSgUC4cBpYDKwxPr6QmCO48MTQghhq1YTutb6JPAicBxLIj8PZAOlWmuTdbV8oIuzghRCCNE6pbW+8gpKdQA+Bm4HSoGPsFyZ/8Za3YJSKhn4wlolc+n284B5AEajMW3x4sV2BVpRUUFkZKRd23ozfyy3P5YZLi53VVVVs+uEh4e7MiSn88f/tT1lnjRpUrbWOr219Wxp5TIFOKq1PgOglPoPMBaIVUoFWq/SuwInm9tYa70AWACQnp6u7R372B/HTQb/LLc/lhkuLre/tHLxx/+1M8tsSx36cWCUUipcKaWAa4A9wGrgFus6c4FPnRKhEEIIm9hSh74ZSxVLDrDTus0C4Gngl0qpQ0A88IYT4xRCCNEKmzoWaa1/Dfz6ksVHgBEOj0gIIYRdpOu/EEL4CEnoQgjhIyShCyGEj5CELoQQPkISuhBC+AhJ6EII4SMkoQshhI+QhC6EED5CEroQQvgISehCCOEjJKELIYSPkIQuhBA+QhK6EEL4CEnoQgjhIyShCyGEj5CELoQQPqLVhK6U6quU2t7kp0wp9ZhSKk4p9bVS6qD1dwdXBCyEEKJ5tkxBt19rPVRrPRRIA6qApcAzwCqtdW9glfW5EEIIN2lrlcs1wGGtdR4wG1hoXb4QmOPIwIQQQrRNWxP6HcD71sdGrfVp6+MCwOiwqIQQQrSZ0lrbtqJSwcApYKDWulApVaq1jm3y+jmt9WX16EqpecA8AKPRmLZ48WK7Aq2oqCAyMtKubb2ZP5bbH8sMF5e7qqqq2XXCw8NdGZLT+eP/2p4yT5o0KVtrnd7aeoFt2Oe1QI7WutD6vFAplaS1Pq2USgKKmttIa70AWACQnp6uMzIy2nDI72VmZmLvtt7MH8vtj2WGi8udnZ3d7DppaWkujMj5/PF/7cwyt6XK5U6+r24BWAbMtT6eC3zqqKCEEEK0nU0JXSkVAUwF/tNk8R+BqUqpg8AU63MhhBBuYlOVi9a6Eoi/ZFkJllYvQgghPID0FBVCCB8hCV0IIXyEJHQhhPARbWm2KIRwM39pzijsI1foQgjhI+QKXQgHkatn4W5yhS6EED5CEroQQvgISehCCOEjJKELIYSPkIQuhBA+QhK6EEL4CEnoQgjhI6QduhBOJu3ThavIFboQQvgISejC+9XXuzsCITyCJHThvUwm+NnPICQEwsKga1eYMAEOHnR3ZEK4ha1T0MUqpZYopfYppfYqpUYrpeKUUl8rpQ5af3dwdrBCXFBejp51A+/87Tw3JW/l1bHvc2LsHbB3L4wfD7t2uTtCIVzO1iv0+cBKrXU/IAXYCzwDrNJa9wZWWZ8L4XynT3N29HXc/uWPuJd3WFeVxiOr5nDVhy8yoVseRToRMjIgJ8fdkQrhUq0mdKVUDDABeANAa12ntS4FZgMLrastBOY4K0ghmsq5/U8M2fM+Sw0383//BwUFsG8f/OEPsGVXOPf02Yw5IgomT4ZDh9wdrhAuo7TWV15BqaHAAmAPlqvzbOBR4KTWOta6jgLONT6/ZPt5wDwAo9GYtnjxYrsCraioIDIy0q5tvZm/lLuqqurCY5PJRGCgpUVteHj4RetFbN3Oj56aQnFEF37758P07Vtx0euffZbESy/15cFbc3lt+ShKhw5l1x/+4PwCcHEZbHFp2Zr+r9u7L2/hL+d3U/aUedKkSdla6/TW1rOlHXogkAr8XGu9WSk1n0uqV7TWWinV7CeD1noBlg8E0tPTdUZGhg2HvFxmZib2buvN/KXcTdtqFxYWYjQagUvaajc08MrdW9nNID55q47Zt15+fk+caLlif2PxEO584E0mLbiTjJoamDHDpWWwxaXt0Jv+r5vbl6qtRdXXYw4Lg4CAK+7LW/jL+d2UM8tsSx16PpCvtd5sfb4ES4IvVEolAVh/FzklQiGsCl/9kF+fmseMlFPccEtws+soBf/8J/TuDXctu52Sq4fDY49BXZ2Lo3WsqO++Y8j06QzLyCBt5EiGjhtH57/9DVr5hi38S6sJXWtdAJxQSvW1LroGS/XLMmCuddlc4FOnRCgEQGUlz/y3gWoVzvwPklCq5VUjI+H996GgQPH3Me/A/v3w17+6LlYHi/viC3o/+ih1SUmceOwxTv34x5SNHk3SW2+R9K9/uTs84UFs7fr/c2CRUioYOAL8CMuHwYdKqfuBPOA254QoBGx6fAlvV8/lmR/k06dv11bXHzbMUsvy16/78uT0Gwj57W/hBz+Ajh1dEK3jGN95h67z51OelsahP/8Zs7XuNf9EEMNCf8dVCxbQEB4OXlrlIhzLpmaLWuvtWut0rfUQrfUcrfU5rXWJ1voarXVvrfUUrfVZZwcr/JTZzP8uTKZTcAnP/r31ZN7ol7+EwkJYnPEPKC+HV191YpCOF/Xdd3SdP59z11zDwVdfpdQcw4cfJnL33f2Zc+MQxu58l8wRPyb5lVfgnXfcHa7wANJTVHi8E4s3sLJmIg9cV0BbGgdMmQKDBsFL7yehZ90A//gH1NQ4L1AHUvX1JL/4IjVdu3L0d7/j9NkIbrxxEC+8cBUADz98kro6A9O3v8b8Lr9HP/kUVFe7OWrhbpLQhcd7+48FmAngR8/3atN2Slmu0nNz4dsJv4HiYnjvPecE6WAd33+fsGPHOPHEE5iDQ3j++W7U1Sneemsfixbt5f77C1i0aC9paeU8dvI5/qfwJ/DGG+4OW7iZJHTh0cxnS3lz13Cu6bqPq/uHtHn7u+4CoxFe+jYFhgyBV17x/JYhJ0+S9PrrlI4fT9m4cXz6aTybNsXwyCMnGTy48sJqHTqYmD//EDOml/CCeppjzy+C2lo3Bi7cTRK68Gjf/nYdx3R3Hng4yK7tQ0Lgpz+FFSsU+2//H9i5E1avdnCUDvbUUyiTiROPP05BQRAvv5xMenoZN9985rJVDQZ45NGTBAQH8FzBT+Htt10fr/AYktCFR3v93VDiAkqZ84ur7d7HvHmWxPd+1WxITLRcpXuq7Gx47z0K77mH2i5d+f3vu2M2w69+lYehhXdrx471PPYLA4v4Adt+u0yGE/ZjktCFx6rens/SsxO4Z8wRQsOu0PC8FUajZVTdj5YGwsMPw/LlnjvE7ksvQVQUBffcQ05OJJs3R/PTn56kS5crd4x6+hlFfHQdT59+FN5910XBCk8jCV14rG//UUEdIdz/u27t3tctt8CePbBn8s8sl+uvv+6ACB0rpKgIPvwQHnwQc2QkS5cmEBlpYs6c4la3jYmB534TxNdM46vnt7ogWuGJJKELz6Q1n+/oT2r0QQZnxLd7dzfeaGn18vHaRLjuOvj3vy0TZLhRdnb2RT8dP/wQbTazMyOD0tIAVq3qwMyZZwkNbf0mbnZ2NiNHbSM5uoTnD9/OrqVL2zy2jPB+ktCFR6rYdpot9ancONEx/dU6d4axY2HJEuC++ywjeK1c6ZB9O4KhqoqrVq7k3OTJ1HXuzPLl8dTXG7jppstvhLYkOFgze3Yxa5lI5cdZToxWeCpJ6MIjbV1cDsCsn7a/uqXRLbdY2qQf6DXTMgTAm286bN/tFf/ZZwRVVlJ0991oDUuXJjJkSAW9erWtI9TUWyx17V9+IROI+SNJ6MIjrc7uTHLgSYZM6+Swfd50k+X3kk+D4J574LPPoMgDBgltaKDj++9zrl8/KgcPZtu2SPLyQrnxxtbrzi/VpUsd6V2O8sHZawk5fMQJwQpPJgldeBxzcSWrz49gSp99VxxVsa2Sk2HUKGu1y49+ZKlDX7TIcQewU8zGjYTm53N0jmXSr//8x3IzdOpU+6qbpt1YzR4GcuqDvY4MU3gBSejC4xxZWkY14YyZaXb4vm+5BbZtg8OhA2HECEu1i5t7jiYsXUp9fDyFo0dz/nzbboY2Z/KcOgKViS++Nrq9bMK1JKELj7Plu0SiVDkD5zi+Hvjmmy2/ly3DcpW+axdkue8GYtCZM8Rs2EDJ9dejAwNZvz6G+noD111XYvc+Y2MbyOh5kA/Lr6chZ4cDoxWeThK68CjmOhPfFg4nIymX4FAH1rdYde8OffvC118Dd9wBoaGwcGFrmzlN/LJlqIYGiq3VLevWxRAfX0///m2bU/RSU+8wcYourHlR2qT7E0nowqPkfZrPaZ3EhInnnXaMadMgMxNqw2Jh9mzL9EbuGNTKbCZh2TLK0tOpTU6mvl6xaVMM48eXttjN31ZjZ9QTZajk3RVxjolVeAWbThul1DGl1E6l1HalVJZ1WZxS6mul1EHrb2knJdpt43KFgQbS7nZeIpo2zTJ0+IYNwNy5cPYsfP65047XkqitWwk5efLC1fmuXXFUVgYwfnz7P8xCQzVT+uzjs7IJmI8db/f+hHdoy3XAJK31UK1141TrzwCrtNa9gVXW50K0S+bB3qRH7CS6k32jK17a+7LpT6OJEyEwEL76Cpg6FTp1svQcdbGETz7BFBND6aRJAGzd2pHgYDMjRpQ7ZP8jp9RTTCI5r+c4ZH/C87Xni91soLHycSEwp/3hCH9WefQ82+sGMrxPnlOPExUFY8ZYE3pgoGWu0c8/hzO298psr4DSUmJXr6Zk5kx0SAhaw+bNHRk+vJywMMe07km7IRSFmZWfeMcsTaL9bE3oGvhKKZWtlJpnXWbUWp+2Pi4AjA6PTviV3UvPojEwYPyVRxZ0hGnTLM0Xi4qwVLuYTJa6dBeJX74cg8lE8ezZABw9GkpBQQTjx5c67Bgd4hpIS8hj5b7u0NDgsP0KzxVo43rjtNYnlVIdga+VUvuavqi11kqpZhu8Wj8A5gEYjUYyMzPtCrSiosLubb2Zv5S7qqqKzeuCCKcS49ggCgsLAdpc9qqqlluHrFix4sLjiIgYYCwvvridmTNLSevdG/72N7KHDLEn/FaPfRGt6btkCef69SMvKgoKC1mxwjLee//+hyksdNwVdVq/cv61fjyrX3kbldbTYft1FH85v5tyZpltSuha65PW30VKqaXACKBQKZWktT6tlEoCmu1DrbVeACwASE9P1xkZGXYFmpmZib3bejN/KXd2djbfnY5mdNwuwiKDMRotX/jS0tLavB9bJCRATIyJvXu78sILQy3TGj32GBmJiTBwYJvjb8uxI7ZvJzI/n2O/+tWFcm7b1pWrrz7PwIExQIxdx2/OD57uyT/XB1CyNYZbHs9w2H4dxV/O76acWeZWq1yUUhFKqajGx8A0YBewDJhrXW0u8KlTIhR+oTS3mP2m3owc7Jp67IAAGDGijM2boy2dKe+6y1Kf/tZbTj924tKlNEREcG7aNABKSwPYuTOCESMcP67MqBmxxASUszIz1OH7Fp7Hljp0I7BeKbUD2AJ8rrVeCfwRmKqUOghMsT4Xwi65n1padgy9tu0TQdtr5MgyzpwJZs8eLFPTzZ5t6WTkxDbpAeXldPjmG85On445LAyALVuiMZsVw4c7/sMsMBCm9Mrjy8IU9Pkyh+9feJZWE7rW+ojWOsX6M1Br/bx1eYnW+hqtdW+t9RSttWMGrhZ+aevWaDqoc3Sf5LruDKNGWRLcV19ZF8ybB8XF8MknTjtm3MqVGGprL7Q9B8jKiiIiooHevZ3TmWrGrCDySWbPOzLhha+TnqLC7bRZs6GgP2M77sYQ4Pju/i3p1KmeLl1qWbvWumDKFMvYAAsWOOeAWpOwdClVffpQ1b//hcVZWVEMHVpBQIBzBtKa/uPuAKx8/5xT9i88hyR04XaHvzzEcXMyI4a5PuGkppazdi2YzVjmGn3wQfj2W6dMIh2+dy/hBw5Q3DgfHnDmTBDHj4eSnu6YzkTNSe4VwsDIPFZuk5bFvk4SunC7b/+dD0DKrHCbt2mtN6itUlMrOHsWdu+2LvjRjyx3TJ0wiXTHDz6gISyMszNmXFiWnR0J4NSEDjBlaDHrq1OpzStw6nGEe0lCF263al0QnQ2n6Twi1uXHTk21JNIL1S5JSTBrlqW1S53jOjgFFhfT4csvKZk1i4aoqAvLs7KiiIw00adP+0ZXbM3E66OoIYysf+9x6nGEe0lCF26lG8xknurDmKT9Dp2dyFadO9fRtSusWdNk4bx5lmEAPnVcS9yOH32Eamig6M47L1qenR3FsGEVBAQ47FDNGn9vDwDWflHp3AMJt5KELtzq4MrDFOmODBvqniZ1SlkG61q7tsnkPtOmWW6OvvyyQ2b8UTU1JC5ZwvkJE6hNTr6wvLAwiBMnnFt/3ighKYiBEUdZs0uG0/VlktCFW61bfBKAQdNtrz93tIkTobAQDhywLggIgCefhE2bLAOnt1P8ihUEnj9P4V13XbQ8K8tS9eKKhA4wYUAxG8qHYCp23HgxwrNIQhdutW6jgQRVQudRrq8/bzRhguX3RdUu991nGVb3f/+3fTs3m+n4/vtU9utHRWrqRS9lZUURHW2id+/q9h3DRhOmh1NBFNvf3eWS4wnXk4Qu3Grd8W6MSzqMMrihAt2qTx8wGpvcGAXL1HRPPGFpwrhpk937jt64kbCjRym66y4uvUmQkxNFamp5u2cnstWEH1oGAFuzzHmzQQn3koQu3ObU1pMcMXVj/Ag3TP/WRGM9+po1l1SZ//jHEBcHzz9v345NJrr89a/UJiVxburUi146fTqYkydDXFbdAtC5Zxi9Qk+wdnu0y44pXEsSunCbdf8+CsD4mzu6ORJLtUt+Phw92mRhZCT84heWyS+2bWvzPhM++YTwQ4fIf+wxdNDFMzA1tj9PS6toT9htNqF3AevODcRc4dxmksI9JKELt1m32kQEFQy7tZe7Q2HiRMvvi+rRAX72M4iOhmeeaVOLl4Dz5+ny979TnpZG6eTJl72+Y0ckkZEmevZ0Tf15o4lTgjhHHLsWSz26L7J1ggshHG7doSRGxx8gMCS19ZWdKDs7G7MZYmJS+OSTUoYMsUyBl5aWBrGx8Ic/WBL7P/4BDz9s0z47L1hAQHk5Jx5//LK6c7Ak9CFDKp1ef35p79mEEfUArF1awpAHnHts4XpyhS7covTIWXbW9mb8UNdWObTEYIAhQyrYvj3y8hd/8hOYPh0ef7xJ28aWhR4+TOKSJRTfeCPVffpc9npZWQBHjoSRkuL6sht7B5EceJK1We5rJiqcRxK6cIsNbx9AY2D89Y6bnae9UlIqOH48lHPnLvniqhS8+aal5csPfgD19S3uI6ioiJ5PPEFDRASnHnqo2XVycyMuHM8dRnY5zNqifuh6k1uOL5xHErpwi3VfVhFEHSPv7evuUC4Yav22sGOHJeFeNPDX6dMcfvpp2LoVfvnLZpN6YHExfR56iKCSEg698gqmDs2P7b5jRyQBAZpBg9zTDT8lpZJCjBxZsa/1lYVXsTmhK6UClFLblFLLrc97KKU2K6UOKaU+UEoFOy9M4WvW7UkgLeoA4XGeMzVa//5VBAWZ2bGjmWoXoHTKFIpuvx3++lcYPhxyciwvmM2wZ48lmZ85w8G//IXKK0w2vWNHJH37VhEa6pzxz1vTf5rlA2vjf0675fjCedpyhf4osLfJ8z8BL2utewHngPsdGZjwXTXnqsmq6Mu4/iXuDuUiISGa/v2rWkzoACeefBKWLrWMFTBiBIwaZblxOnAgwYWFHPrLX6gcOrTF7evrFbt3R7itugUgOT2caFXGxo3u68wlnMOmhK6U6gpcB7xufa6AycAS6yoLgTnNby3ExbLf208dIYyd4nk35lJSKtizJ5yamiskuzlzYM8ey2QYQUFw773w+uvsWbyYimHDrrj//fvDqK01uDWhBwQqRnU8wsa8Lm6LQTiHrVforwBPAWbr83igVGvdeFclH5CzQ9hkw3LL9LNj7unp5kguN3RoBSaTgb17I668YocO8Pe/w7p1liqY+++nrkvrb4HGVjRD3dy6Z8zQanbW96XsgEx44UtabYeulLoeKNJaZyulMtp6AKXUPGAegNFoJNPO0esqKirs3tab+WK5V28OpFfgEfYUHGePNZ9UVX3fc9FkMlFYWAjQYtmbru9ISUlngV6sX6/p3Lmw2XXaE9OWLV0wGqtoaDhJ4SW7b1puZ0vsbUJ/aWDx89/Q50ddXXLM5vji+d0aZ5bZlo5FY4EblFIzgVAgGpgPxCqlAq1X6V2Bk81trLVeACwASE9P1xkZGXYFmpmZib3bejNfK7duMHNr6Vmu77X/onI17QBTWFiI0WiZ/zItLa3Z/dgz3ZwtjEbo3r2aw4c7YTQ23wrF3pi0hv374xkxouxC+ZpqWm5nG/c/g/j5Xxs4fTCReW48v3zt/LaFM8vcakLXWv8X8F8A1iv0J7TWdyulPgJuARYDcwHHTe8ifNaBL49SrHsydqxtTeaclbivJCWlktWrYzGbcWhPzvz8EEpKgtxaf94oOjGEwREH2bjXfcMWC8drz+n6NPBLpdQhLHXqbzgmJOHLNnxo+SI39lbPveUydGgFZWWBHDvm2CaVje3b3V1/3mhMrzN8V9qPhir3jnYpHKdNCV1rnam1vt76+IjWeoTWupfW+lattZwVolUbNkKcOkvf6d3dHUqLGq+gmx0GoB1ycyOJiGjg6qtrHLpfe43JCKaMGHZ/LB2MfIUMziVcasOxrozpeAhDwAh3h9Ki5ORa4uLq2bEjkptuKr7sdXurgXbujGDQIOcPyGWrMXd2g/mw8dMzDLnH3dEIR/CQU0v4g+K9Z9hffzVjh7l2yNi2UsoyUFdjFYkjVFQYOHw4jCFDPKO6BaDHiESMAWfYuDWo9ZWFV5CELlxm4zuHARg7y/Nnnk9JqSQ/P5SSEsd8id29OwKzWTFkiHvGb2mOUjCm8zE2nryqTWO9C88lCV24zIZvqgmijvQ7e7s7lFY11qNfaRiAtti503K1764BuVoyJq2Oww09KMzOd3cowgEkoQuX2bAvjrTIA4R18JwBuVrSr18VwcEtD9TVVrm5kVx9dTVRUQ0O2Z+jjLkhAYDvFh9zbyDCISShC5eoPV9DVnlfxvbzrAG5WhIcrBkwoNIh9ehmM+zaFeFR9eeNUm/tSRB1bMyURmq+QBK6cIns9/ZTSyjjpnj+1XmjlJRK9u1rZaAuG+TlhVJWFsjgwZ5V3QIQGhlIatQhNu2Pd3cowgGk2aJwifWfWQbkikmtcUvvT3ukpFSwcGEn9uyJIDXV/qvrxhmKPOmGaFOj+57ln1mp1JdVExQd5u5wRDvIFbpwiQ05YfQJOkr01Y7trONMjVUk7a1H37kzgqgoE926eUaHokuNmRxCNeHs+HC/u0MR7SQJXTidbjCzoag3Y7s3O36bx4qNbaB79+p216Pn5kYyeLDndCi61Og7ewCw8TPvuL8hWuahp5jwJftXHqVExzNunLsjabuUlEpycyMxm1tftznl5QEcOeJZHYou1XVoAl0DTrMpR2aR9HaS0IXTrf/AOiDXbZ47IFdLUlLaN1DXrl2Wq3tPvCHa1JiueWw61V06GHk5uSkqnKLpjc9Vq+tIUMWUx5eg8K55LJt2MLJnUK3c3AgMBs3AgZ6d0Eenm/gwL5lVC78idvDlLV5aGgNeeBa5QhdOt7WgJyPi9qEM3pXMAa66qpYOHertHnkxNzeSXr2qiYy0s87GRUbfkAjA/q/K3ByJaA9J6MKpzh8q57CpB8MGeOcNN6UsV+n2JPSGBs/tUHSpYbf0JIQaduV4Tz8BcTlJ6MKp9q8oBWDgBO+t3Rs6tIKTJ0M4c6ZtoxIePhxGZWUAKSmeXd0CEBweSHr0QXKOe999DvG9VhO6UipUKbVFKbVDKbVbKfVb6/IeSqnNSqlDSqkPlFJyi1xcZseWYEKppse0BHeHYrdhwyxX2Nu2te0q/fsORZ5/hQ4wun8p26oHUF8mwwB4K1uu0GuByVrrFGAoMEMpNQr4E/Cy1roXcA6433lhCm+Vk9eF1Mi9BEV47xV6375VhIU1tDmh79gRSXx8PZ071zkpMscac00YdYRwdOUZd4ci7NRqQtcWjZcYQdYfDUwGlliXLwTmOCVC4bVqz9WyrXoAqT1OuzuUdgkMtDQ7bGs9em5uJCkpFSgvuRc85p6eAOxe59k3cEXLbKpDV0oFKKW2A0XA18BhoFRrbbKukg9I5Zu4yMHlxZgIImWsZw0Za49hwyo4dCiM8vIAm9YvLg7k5MkQr6luATD260DPwGNsk4G6vJZN34O11g3AUKVULLAU6GfrAZRS84B5AEajkczMTDvChIqKCru39WbeWu6qqiqyVtVjoIH40WYKCwtt3tZkMrVpfVfo1q0erTuTmVnPiBGtV0ls2NAJgOTkExQWltp0DFeWu6VzKrVjHatPpVJw6jtUgKHV9dvLW8/v9nBmmdtUsam1LlVKrQZGA7FKqUDrVXpXoNmBOrTWC4AFAOnp6TojI8OuQDMzM7F3W2/mreXOzs5mxzFNSuherhqY3KZtCwsLMRqNTorRMP0mAAAcX0lEQVTMPhMmKAIDzRw7lsysWa1/sc3L60JwsJkxY0IJCrKtLK4sd0sdhb5LW8ZHpxLQhwPoNL5jq+u3l7ee3+3hzDLb0sol0XpljlIqDJgK7AVWA7dYV5sLfOqUCIVXqq+oZ2vFQEZ0O+7uUBwiNFTTv3+VzTdGc3MjGDCgkqAg7+pKP2BqOAB7v/HsibxF82ypQ08CViulcoGtwNda6+XA08AvlVKHgHjgDeeFKbzNkS/OUEMYKaPr3R2KwwwdWsGePa1PeFFTo9i7N9xjxz+/kqSRscSpsw6bek+4VqtVLlrrXGBYM8uPACOcEZTwfrmZlivTfrM6uDkSxxk2rIJ33unE7t0RpKW1fLNz375wTCbDhXFgvIkhQDE8bh9bC64Gzrs7HNFG0lNUOEXWvo4MCDpAdLf2z8npKRoTdGvVLo1Xt954hQ4wrG8JB0y9KDskCd3bSEIXDmeqMbH5/EBGXnXU3aE4VExMAz17VrNtW9QV19u+PZKrrqqhQwfTFdfzVAMnWJpmNg7bILyHJHThcNs+OEAFUQwd7ns31tLSytmxI4L6+ubr0RsaLFfwaWnlLo7McXpOjyOYWnZuldE8vI0kdOFwa5cUAdB/VqybI3G89PRyamoC2L07vNnXDx8Oo6Ii8ML4L94oOCqIoeH7yDkmfQW9jSR04XBrt4bRM/AocX2vXDXhjVJTy1FKk5XVfNlyciIvrOfN0nqcIqd6AHXnZaAub+K9IyYJj2Q2mVlX1Idrk7MB32nh0ig2toHevavJyorigQcKLnt927YoOnWqpVMn726uOWRkPXW7Qzi8opj+d3a5aAaqS8lsRp5DrtCFQ+V+fJBzuoPXX6FeSXp6Obm5kdTWXlyPrrWl/tybq1sa9Z8dh8LMjkzv6hjl7yShC4da9a5lZMUhN8W4ORLnSU8vp67OwM6dFzfJzMsL4ezZIFJTvT+hR3YJY1DIfrYeSHJ3KKINJKELh/rmuwj6BR8mboDvJvTU1HIMBk129sX16I3D6w4b5hvfTkZ1z2Nz+WDqy71jPHchCV04UF1FHWuLB3BN32bHafMZkZFm+vWrYuvWixN6Tk4UcXH1dOvmGzcSh42po4YwjqwocncowkaS0IXDfPfWXqqIYMp1Ie4OxenS08vZtSvionFdGuvPvWVCi9b0nx0HwPbVbg5E2EwSunCYVR+dxUADGQ/ZPFy+10pPL8dkMlyoZjl9OpjTp0N8proFIKprGINC9rFlv9SjewtJ6MJhvtkeT3rEXmK7+W79eaOhQysICPi+Hr1xfBdfuCHa1Khux9hSPghThdSjewNphy7s1rRtclVhNZvLR/GTISvJzvaNOuQrCQ83M3Bg5YV69JycSCIjTfTs6VvDHQwbXcvrB8I5suIMfW6TnqOeTq7QhUPsXVpCA4GkTfKfCYZHjixjz54Izp0LYMuWaNLTywmwbcpRr/F9Pbq0R/cGktCFQ2SvCyKUavrMTnR3KC4zZkwZZrNi5co4Tp0KYcQI36k/bxR9VTgDg/ezZX8nd4cibCAJXTjEhqM9GRWdS3C0/4zQN2BAJdHRJr76ynIVO3JkmZsjco5R3Y+ypWwQpkrvHs7AH9gyp2iyUmq1UmqPUmq3UupR6/I4pdTXSqmD1t++N3CHsMn5Q+XsruvLqIGumbHeUwQEWJL4/v3hGI11XHWVb947GDa6lioiOLxc2qN7Oluu0E3A41rrAcAo4KdKqQHAM8AqrXVvYJX1ufBDOYstEyGkzfS/e+yjRpVRV2egX79Kn2l/fqnBN8djoIGtX7s7EtGaVhO61vq01jrH+rgc2At0AWYDC62rLQTmOCtI4dnWbexAkqGA7tM7ujsUl0tMtDTnCwnx3ZvBEZ3DSA/fybp9PdwdimhFmy6plFLdsUwYvRkwaq1PW18qAIwtbDMPmAdgNBrJzMy0K9CKigq7t/VmnlzuqqoqTFUmMosmMrPLRorOOKYlhMlkorDQO6pvsrMt7c+PH1ftjtmV5W7pnKqqqmp2+cheNfw1dw7HspYSlhxm075s4cnnt7M4s8w2J3SlVCTwMfCY1rpMNfl+qbXWSqlm381a6wXAAoD09HSdkZFhV6CZmZnYu6038+RyZ2dns+udk5QRw4SpdRiNXR2y38LCQozGZq8PPM6ePUl06FDPwYMdiIxMIiLC/it1V5a7pTHMWxr3fMwNBbyaayD/WwPjnro4xvaMh+7J57ezOLPMNrVyUUoFYUnmi7TW/7EuLlRKJVlfTwLkjokf2vRlIMHUMuSOBHeH4nI1NYrt2y3zhzY0qMsG6/IlV89MpIM6x+aNke4ORVyBLa1cFPAGsFdr/VKTl5YBc62P5wKfOj484elWH+rLuJhthCaEujsUl9u+PZL6egMzZ54lPLyBTZt8d8iDgOAAJnbcwZqTg9Bm6WTkqWy5Qh8L3ANMVkptt/7MBP4ITFVKHQSmWJ8LP1K45SwHTD2ZkHb5VGz+YPPmaIKCzAwfXs7w4eVs2BCN9uFcN3pEKad0EvnfypdxT9VqHbrWej3QUoOsaxwbjvAmWz6qBGD4rRGtrOmbNmyIISWlgrAwMxMmlLJmTSz794fRr59vjefSaNitMfAZZH1WR/IUd0cjmiM9RYXd1mYn0TfoEMbhce4OxeXy84M5ciSMCRPOAzBxYikGg+bbb323f13cgGgGBB1g/c7O7g5FtEASurBL+aly1pcNY1KvA+4OxS3Wro0FYMIES6eq2NgGUlPLWb061p1hOd34XofZVJZC7Vnf7BXr7SShC7usfHEX9QQzeobJ3aG4xZo1sfTsWU3Xrt+PEz5pUilHj4Zx9Kjv3iAeeU09dYSw64Pv69Gzs7Ob/RGuJwld2OXDD6GjKqLvLf43Ct/58wFs3x554eq8UUaG5bkvX6X3v9VINOdZ85V/3jfxdP43+IZos0uvtmqKa/j85DDu6LmOwBD/a3++cWMMDQ2KiRMvTuhGYz2DBlWwenUs993nmy1/giICmdY5my/zU3mk5hCBoZJCPIlcoYs22/pGEdWEM/lG/5yWbO3aGOLj6xkw4PJu8pMnl7J3bwSnT/vuMMKTppRTouPZ95F3DM/gTyShizb7+ttEOhtO+2V1S12dYuPGGCZMKMXQzLtn0iTfr3ZJubcjYVSxdrnvfmh5K0nook0qC6r5pmQ41/XajiHQ/06fnJxIKisDLqtuaZScXEvv3lV8+63vJvTQ2GCuSczii6NDMZt8d5RJb+R/70jRLllvnqGOECbd4p9v5DVrYgkNbSA9veXp5iZPLmXHjkgKCoJcGJlrTZ54ltPmThz61DfvFXgrSeiiTb5cbSTZkE+vG/xv7HOTCb79tgOjR5cRGtpyH/+ZM0vQWvH55/EujM610n6YQBB1rPtEUognkf+GsFnlySpWn0vnur65flndsnVrNCUlQVx77dkrrtelSx1paeUsXx7vs2O7RHQKI6NDNl8ckMG6PIn/vSuF3Tb8o4R6gpl0m3++gVesiCMqysS4cedbXXfWrGJOnAhlxw7fba89eUwhxxq6kfeVDNblKSShC5t9+G1vBgbv4+rrvGPyCUeqqjKwenUsU6acIzi49Q+0yZNLCQtr4LPPfLed/sj74giknm8WuTsS0UgSurDJ0RWFbK8dyG3j96EMPjob8hVkZsZSUxPAzJlXrm5pFB5uZsqUc3zzTQeqq33zbRbdLYIZCd/x8b7hmKr9cwgIT+ObZ5pwuGVvBhJGFZMe890rzitZsSKOzp1rSUmpsHmbWbNKqKwM8Ok26bNuPMsZnUjOv063vrJwOknoolVVhdX859gY5iRvJDwprPUNfExxcSBbtkQzY8bZZjsTtWTo0Aq6dKll+XLfbe2SMrcTSYYCPv0s0d2hCGybgu5NpVSRUmpXk2VxSqmvlVIHrb99dxBowbpXi6kkkht+6J9Dpq5cGYfZrJg5s6RN2xkMcP31JWzZEk1eXoiTonOvwNBAbhmUxdfnRnFub+s3i4Vz2XK98TYw45JlzwCrtNa9gVXW58IHabPmw9W9GRS8j16z/O9mqNbw+efxDBhQSffubf9Au+mmMwQFmXnvPd/92814KAQzAaz6m+3VUcI5Wk3oWuu1wKV3gmYDC62PFwJzHByX8BBZ7+xlR+1Abp/onzdDs7MjOXgwnJtuOmPX9vHxJq67roTly+M5d843RyY0johjbGQ2H2UNljbpbmZvHbpRa914F6QA8N3LDz/3/547TxRlTHjUP2+GLlpkJDa2nhkzbGvd0py77y6ittbAkiW+W89845TjHDb1YN8Hp9wdil9r9yWD1lorpVr8WFZKzQPmARiNRjIzM+06TkVFhd3bejN3lrtwUylL8m/gof6fUKFCqSgsc8lxTSYThYXuH5r15Mlw1q+P4fbbD3H+fAHn7awiDg+H4cMTWbw4nunTcwkJaX4cHFeWu6Vzqqrq8iGBbdHzVgMdPjnLojdj+Pnkwiseoyl/fF87s8z2JvRCpVSS1vq0UioJaLGrmNZ6AbAAID09XWdkZNh1wMzMTOzd1pu5s9z33LeeMKq59X8TiDa6rsdjYWEhRqP7v/S9/XYygYGaH/6wioSE9sVz//3neOihjmRn9+fGG4ubXceV5U5LS2t2ud1Txxnhh6kbeDlnFg/vW0mXiYktHqMpf3xfO7PM9la5LAPmWh/PBT51TDjCUxxalcd7R0fzcNpWorv5bvf1lpSVBfDZZ/FMn36WhIT2d5pJS6ugf/9K3n3XiNlHB6q84dlowqlk0UuWUSZlrlHXs6XZ4vvAJqCvUipfKXU/8EdgqlLqIDDF+lz4kP97+DjB1PHEmwPcHYpbLF2aQE1NAHfd5ZhxSpSCe+8tJC8vlC++iHPIPj1NVLdI7um7ho9PTuRMjv33HIT9bGnlcqfWOklrHaS17qq1fkNrXaK1vkZr3VtrPUVrLf89H3JsfT7/PjiKB4dsodMQ/xsmt6ZG8cEHHUlPL6NPn2qH7feaa84xYEAlr73WhZoa32wxdPN/h2LAzId/8tGvIR5OeoqKy/zq3mMYMPPU633cHYpbvPeekaKiYB580LHd2Q0G+MUv8iksDGbRIvffI3CGuIEx3N5tDe8dnkjpAdfcRBffk4QuLpL5ynbePTqOp8dvouvwJHeH43IlJYG8/XYnJk4sJS3N8R1lhg2rYNKkcyxc2IniYt9sl37Hk5p6glj8mxp3h+J3JKGLC+oq6vjJ01H0CDzOf30y0t3huMW//pVEba2BRx7Jd9oxfv7zk9TWGvjnPzs77Rju1GlUPD/osYq3Dkzl+DcyVrorSUIXF7x860b21vXkr88VEhbnf4NwHT0aytKlidx88xm6dXPeuDVXXVXLbbcV8emnCezZE+6047jTfS+FEqvO8/9+b5Teoy4kCV0AkLchn9+tHM6cpO+Y+evh7g7H5bSG+fO7EBpqZt485/d2fPDB0yQm1vPccz18crz0yOQInrxuHZsqU1nzhxPuDsdv+N6ZJNrMVGPiR9dbxip55aOubo7GPT7/PI7162N58MFTxMY2OP140dEN/Pa3RzlxIoSXX/bNv/nk57owInw7L3w6lqrTjmstJFomCV3w7MT1rC4dxmsPbKPbWN9MLldy/HgIf/rTVaSmlnPnna6r801Pr+Ceewr5z38SWbMmxmXHdRVDoIGn/vsUxTqe135SI1UvLiAJ3c8tffo7XtiSwUMD1jL3X+PcHY7L1dcrnn22B0FBmt/97igBAa49/sMPn6Jv3yp+//tuFBeHuvbgLtB9hpFHhq1g8YkpfP0/UvXibJLQ/diBL48y94UBjIjYxSub/LNVy2uvdWbv3gh+9as8OnWqd/nxg4I0zz9/hPp6A889N5zSUhd/orjAnX8zMjV2E79eeR0HPylwdzg+TRK6nzqSeZxp1wURrOr56OsOhET75ow6V/LJJ/G8804nbrrpDJMmlbotju7da3nppUMUFITz2GO9qKryrbdlQHAAzy000TXgFI//IZWCXGnK6Cy+deYImxxalcfEKYGUmSP46p0irhrdBfCvwZRWrIjj+ee7MWbMeZ54wv1VAWlpFTz99Db27IngySd7UlfnW0MDRHQJ58/P7+WcOYaZo0o4s7f5ESdF+0hC9zP7VhxhwrQQqs0hrF5cROrd/d0dkst9800sv/lNd9LSynnhhcMEB3vGzbrRo4t47rk8Nm+O5qGH+lBS4ls9Sa+a0pFXH/qafdXdmDD0PPlbHTu0gpCE7lcWP7KRkdfFY9IBrP6ohJTb+ro7JJfSGhYt6sizz17N4MGVvPTSYUJDPSOZN7rhhhL+7/+OsH9/OPfe2599+3yrg9fQB7ry5asHOVWXwLjRDZzJkomlHUkSuh+oLKrkvj7ruPPVMQyMPM6WtbUMvtm/Bt4qLw/gqaeu5uWXkxk37jzz5x8kPNwzRwScOvUcb7yxD4D77+/H4sWJmNo/JLvHGP+zFFa/e4pKcygPPJnBW/etkyaNDiIJ3Yc11DXw73nrGdj5HG8fHMuzYzNZe6Y/3cf5T1tzrSEzM4Yf/KA/a9fG8otfnODFFw8TGemZybxRv37VvPPOXtLSynnxxau4557+bNsW6e6wHCb17v5syaxmcMRB7ntrPDckbeVUjrSAaS9J6D6otqyWD3+xkSGRh5n7r3HEBZXy7yc+58b5UezYvcPd4bmE1pCdHcl99/XliSd6ERioef31/dx9dxHKS+43xsWZmD//EC+8cJjy8gAefLAvjz7ai02bon1i1qMeE5L53SdlvHLjGr4pGkzvtCgeG7aGE5tloml7teuui1JqBjAfCABe11rLzEVu0lDXQNY7e1n0lxIW7RzMWT2G3oGHefXuTxj1aFeUwTdH9rtUcXEgK1fGsXx5PIcOhdOxYx3PPXeM668vIdAL7zEqBZMnlzJ6dBmLFnXko4868vOf96Zbtxquu66ECRNK6dmzxms+pC5lCDTw6H8mct2qPJ5/6Dh/2z6G10Zpbuu+gdvuCmTakymExvpehytnsfsUV0oFAH8DpgL5wFal1DKt9R5HBeePLm0mWFVVdWFZ00l3i/eXsPPz4+xYV0bmd6FkFvTjPIMIoYY5yTn8cF4wcdeYCQhOdmn8rqQ1FBUFsX9/ODk5keTkRLFvXzhms2LQoAqeeSaP668v8bgbn/YICzPzwAMFzJ1byDffdOCjjxJ57bUuvPZaFzp3riUtrZyBAyvp37+Knj2rva7Mva7pxlsHu/GbDfm8+PBh3t2VwqI/xBL5h3Kmd9nG6GG1jLw2jtTbehGe4JsjVDpCe65ZRgCHtNZHAJRSi4HZgCT0Fmizxmwy01DXQENdA6YaE/XVlp/a8jpqyuo4vr2ImvMNVJeZqSo1c+ZkNbXVJkpLA/ljxSaOl0ZxrMpIkU4E4gHoHniCW/rkcs20AGb8cgAdeowB2jGDu4uYzZafhgaFyaSor7f8LigIo7IylJoaA1VVBsrLAykrC6C0NJDi4iDOnAni1KkQ8vJCqay09KwMCjIzeHAlDzxwmmnTztK9u/OGv3WnoCDNtdee5dprz1JcHMj69TGsWxfLunUxfPZZwoX1EhPrSE6uxWisIz6+nrg4E7GxJiIjG4iMbKChAcLCvv8JCoLgYKiuNmAwaAICNAaDZZYlV+o2tiuv5nblzxV1rJ6fxcfvVvHVwav5eHlXWA7qp2a6BpykZ1QRVydW0CmxgY5GRccuQcQkBBGdGEJUQghhMcGERlt+gsKDCAoLJCg8iIDgAAKCA1AGL/1K0wqltX2f5EqpW4AZWusHrM/vAUZqrX/W0jbp6ek6KyurzcfqGlbCyRpfmFi3PSeRRqExoDEo65stQGEIVKgW3nVmc/OjBhoMzXcvb2n9i6Jo4XRpulxrddHy73/UhWVms31/i4iIBhIS6jEa6+jevYYePWro2bOaAQMqve6q9FKFhYUYjfZNTac1FBQEs3t3OMeOhZKfH0J+fghFRcGUlARRW2t/ZjYY9IUqncbHTat4lNJNHre+v4AmA+Y0NJgICGj9utJsNmOu15jN2nIhoJX1HaGw733l+nNl5ev5TL8/mczMTDIyMtq0rVIqW2ud3tp6Tq9VVErNA+ZZn1YopfbbuasEwB+7l10otwYagIYLD4A69wXmRC3+rysrLT95ebBli4ujcj6PPMeb3oBtcPzIwh5ZZmeY8QDwAGBfmbvZslJ7EvpJoGkFbVfrsotorRcAC9pxHACUUlm2fEL5Gn8stz+WGfyz3FJmx2pPDdlWoLdSqodSKhi4A1jmmLCEEEK0ld1X6Fprk1LqZ8CXWJotvqm13u2wyIQQQrRJu+rQtdYrgBUOiqU17a628VL+WG5/LDP4Z7mlzA5kdysXIYQQnkW6/gshhI/wmISulIpVSi1RSu1TSu1VSo1WSsUppb5WSh20/u5gXVcppf6ilDqklMpVSqW6O357tVDuW5VSu5VSZqVU+iXr/5e13PuVUtPdFXd7tFDm/2d9nquUWqqUim2yvteXGVos9++tZd6ulPpKKdXZuq5PnOPNlbnJa48rpbRSKsH63CfKDC3+r3+jlDpp/V9vV0rNbLK+Y85xrbVH/AALgQesj4OBWOAF4BnrsmeAP1kfzwS+wNKjYBSw2d3xO7jc/YG+QCaQ3mTdAcAOIAToARwGAtxdBgeVeRoQaF32pyb/a58o8xXKHd3k9UeAf1gf+8Q53lyZrY+TsTSoyAMSfKnMV/hf/wZ4opl1HXaOe8QVulIqBpgAvAGgta7TWpdiGUpgoXW1hcAc6+PZwL+1xXdArFIqycVht1tL5dZa79VaN9cBazawWGtdq7U+ChzCMgSD17hCmb/SWjeO+v0dln4N4ANlhiuWu6zJahF834XR68/xK7yvAV4GnuLiLpteX2ZotdzNcdg57hEJHcun0hngLaXUNqXU60qpCMCotW6cp6oAaOwX3QVoOhFkvnWZt2mp3C3xhXLbUub7sFypgW+UGa5QbqXU80qpE8DdwP9Y1/eFcjdbZqXUbOCk1vrSsZx9ocxw5XP8Z9bqpDcbq5BxYLk9JaEHAqnA37XWw4BKLFUsF2jLdxNfa5LTarl90BXLrJR6FjABi9wTntO0WG6t9bNa62QsZW5xLCQv1FyZfwP8N99/cPmilv7Xfwd6AkOB08CfHX1gT0no+UC+1nqz9fkSLH+QwsavXNbfRdbXbRp2wAu0VO6W+EK5WyyzUuqHwPXA3dYPcPCNMoNt/+tFwM3Wx75Q7pbK3APYoZQ6hqVcOUqpTvhGmaGFcmutC7XWDVprM/Avvq9WcVi5PSKha60LgBNKqcZZi6/BMgzvMmCuddlc4FPr42XAvda74qOA802qZrzGFcrdkmXAHUqpEKVUD6A34FVDVLVUZmWZLOUp4AatdVWTTby+zHDFcvdustpsYJ/1sdef4y2UOUdr3VFr3V1r3R1L8ku1ruv1ZYYr/q+b3g+4Edhlfey4c9zdd4Ob3OkdCmQBucAnQAcsA36vAg4C3wBx1nUVlsk1DgM7adISxNt+Wij3jVhO9FqgEPiyyfrPWsu9H7jW3fE7sMyHsNQjbrf+/MOXynyFcn9sfWPnAp8BXazr+sQ53lyZL3n9GN+3cvGJMl/hf/2OtVy5WJJ4UpP1HXKOS09RIYTwER5R5SKEEKL9JKELIYSPkIQuhBA+QhK6EEL4CEnoQgjhIyShCyGEj5CELoQQPkISuhBC+Ij/D4um3/73cZBpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw(params) ###variationalinference1graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import digamma  ###variationalinference1responsibility\n",
    "\n",
    "def responsibility(z, K, ps):\n",
    "    tau_sum = sum([ps[\"tau\"][k] for k in range(K)])\n",
    "    r = {}\n",
    "    for k in range(K):\n",
    "        log_rho = (digamma(ps[\"alpha\"][k]) - math.log(ps[\"beta\"][k]))/2 \\\n",
    "                            - (1/ps[\"zeta\"][k] + ((ps[\"mu_avg\"][k] - z)**2)*ps[\"alpha\"][k]/ps[\"beta\"][k])/2 \\\n",
    "                            + digamma(ps[\"tau\"][k]) - digamma(tau_sum)\n",
    "                \n",
    "        r[k] = math.exp(log_rho)\n",
    "       \n",
    "    w = sum([ r[k] for k in range(K) ]) #正規化\n",
    "    for k in range(K): r[k] /= w\n",
    "    \n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>610</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>3.969494e-08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>613</td>\n",
       "      <td>0.999999</td>\n",
       "      <td>1.058971e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>613</td>\n",
       "      <td>0.999999</td>\n",
       "      <td>1.058971e-06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     z         0             1\n",
       "0  610  1.000000  3.969494e-08\n",
       "1  613  0.999999  1.058971e-06\n",
       "2  613  0.999999  1.058971e-06"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>500</th>\n",
       "      <td>628</td>\n",
       "      <td>0.315243</td>\n",
       "      <td>0.684757</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>501</th>\n",
       "      <td>628</td>\n",
       "      <td>0.315243</td>\n",
       "      <td>0.684757</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>502</th>\n",
       "      <td>628</td>\n",
       "      <td>0.315243</td>\n",
       "      <td>0.684757</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       z         0         1\n",
       "500  628  0.315243  0.684757\n",
       "501  628  0.315243  0.684757\n",
       "502  628  0.315243  0.684757"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>640</td>\n",
       "      <td>0.000039</td>\n",
       "      <td>0.999961</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>640</td>\n",
       "      <td>0.000039</td>\n",
       "      <td>0.999961</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>640</td>\n",
       "      <td>0.000039</td>\n",
       "      <td>0.999961</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       z         0         1\n",
       "996  640  0.000039  0.999961\n",
       "997  640  0.000039  0.999961\n",
       "998  640  0.000039  0.999961"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rs = [responsibility(d[\"z\"], K, params) for _, d in data.iterrows() ] ###variationalinference1calcr\n",
    "\n",
    "for k in range(K):\n",
    "    data[k] = [rs[i][k] for i,_ in data.iterrows()]\n",
    "    \n",
    "display(data[0:3], data[len(data)//2:len(data)//2+3], data[-4:-1]) #データの先頭、中盤、後ろを表示"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
